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Open quantum system parameters from molecular dynamics 
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We extract the site energies and spectral densities of the Fenna-Matthews-Olson (FMO) pigment 
protein complex of green sulphur bacteria from simulations of molecular dynamics combined with 
energy gap calculations. Comparing four different combinations of methods, we investigate the 
origin of quantitative differences regarding site energies and spectral densities obtained previously 
in the literature. We hnd that different forcehelds for molecular dynamics and varying local energy 
minima found by the structure relaxation yield significantly different results. Nevertheless, a picture 
averaged over these variations is in good agreement with experiments and some other theory results. 
Throughout, we discuss how vibrations external- or internal to the pigment molecules enter the 
extracted quantities differently and can be distinguished. Our results offer some guidance to set 
up more computationally intensive calculations for a precise determination of spectral densities in 
the future. These are required to determine absorption spectra as well as transport properties of 
light-harvesting complexes. 


I. INTRODUCTION 

Understanding the influence of the protein environ¬ 
ment and exciton phonon coupling on the electronic ex¬ 
citation transfer in natural light harvesting systems is 
one of the great challenges of photosynthesis research. 
Unfortunately, a full quantum mechanical treatment of 
even small photosynthetic complexes seems out of reach 
in the near future. Instead, one seeks approaches that 
keep as much microscopic detail as possible and are still 
tractable with reasonable computational effort. One vi¬ 
able approach is to treat nuclear motion classically using 
molecular dynamics (MD) and excitation transfer quan¬ 
tum mechanically, using methods for open quantum sys¬ 
tems. The latter require input parameters that are ex¬ 
tracted from the MD simulation. 

In this article we compare different implementations 
of MD simulations in combination with different ways of 
extracting open quantum system parameters from these, 
to understand large deviations among parameters found 
in this manner in the literature HHS]. As model system 
we consider one subunit of the trimeric Fenna-Matthews- 
Olson (FMO) complex of C. tepidum, shown in FIG. 
since it is most widely studied and with only eight Bac- 
teriochlorophyll (BChl) pigment molecules per subunit 
quite small and tractable. 

To investigate excitation energy transfer or optical 
spectra, knowledge of excited state dynamics is required. 
One has to consider the presence of at least one elec¬ 
tronic excitation in the photosynthetic complex, which 
can migrate among its pigment sites and even be shared 
coherently between several of them. These processes in 
turn also react back on the internal degrees of freedom of 
the molecules and the protein, since a change in electronic 
state induces a redistribution of charge-density and leads 
to vibrations. For these reasons, excited states pose se- 
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FIG. 1: (a) Sketch of the FMO monomer including the protein 
(brown stripes), (b) FMO monomer without protein, showing 
the numbering of the BChl molecules, (c) Structure of an 
individual BChl molecule. 


rious problems for MD simulations. In contrast, ground 
state MD is capable of simulating the thermal motion 
of systems with tens of thousands of atoms with reason¬ 
able resources, and can therefore be used for quite large 
light harvesting complexes. Fortunately, one can obtain 
some of the essential parameters determining the elec¬ 
tronic excited state dynamics also from MD simulations 
in the ground state. 

Typically, electronic excited state dynamics is tack¬ 
led using open quantum system methods, where the rel¬ 
evant electronic levels of the BChls span the system of 
interest and the remaining degrees of freedom form the 
quantum environment. These methods require two es¬ 
sential inputs: (i) A system Hamiltonian, which contains 
the so-called site energies of the BChls on the diagonal 
and the off-diagonal elements of which describe exciton 
transfer couplings between the BChls. For simplicity one 
neglects higher- or multiply excited electronic states of 
the BChls, thus reducing the problem to the so-called 
one-exciton manifold, (ii) Bath spectral densities, which 
describe the couplings of the BChl electronic excitations 
to baths of harmonic oscillators. The question is then 
how to obtain these quantities, which require knowledge 
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of the electronically excited state, from the molecular dy¬ 
namics simulations. 

The method discussed here, is to calculate the (time- 
dependent) electronic transition energy of the individual 
BChls, which are often referred to as ’energy gap func¬ 
tions’, along a ground state MD trajectory. The transi¬ 
tion energies of the BChls depend both on the geometry 
of the BChl and also on the conformation of the sur¬ 
rounding protein, which creates an electric field at the 
BChl locations. The average of the transition energy of 
each BChl along the trajectory is then taken as the site 
energy of the respective BChl. Similarly, the (time de¬ 
pendent) transition dipole-dipole couplings between pairs 
of BChls are calculated along the MD trajectory. The 
main difficulty is to obtain the spectral densities describ¬ 
ing the influence of the quantum environment on the sys¬ 
tem Hamiltonian. Here, as in most previous works, we 
focus on the spectral densities of the transition energies. 
To this end, a classical ’gap correlation’ function is cal¬ 
culated from the transition energies of each BChl. Then, 
using semi-classical a posteriori corrections, the quantum 
correlation function of the environment is constructed, 
see Refs. [SIE]. 

This method has been applied in recent years to various 
light harvesting complexes [IHIl]. The electronic system 
Hamiltonians extracted from these calculations typically 
agree reasonably well with that extracted from fitting 
simple exciton models to experimental spectra. However, 
they show a large spread in parameters, even the ener¬ 
getic ordering of the site energies is not agreed upon by all 
calculations. Furthermore, there are large deviations be¬ 
tween the gap-energy correlation functions and resulting 
spectral densities obtained in different works and their 
agreement with spectral densities extracted experimen¬ 
tally from fluorescence line narrowing techniques [TSHIS] 
is at best qualitative and often very poor. Also funda¬ 
mental questions about the applicability of the approach, 
such as whether a harmonic model is appropriate, remain 
open and require more detailed investigations. 

A crucial element of the scheme is the calculation of 
transition energies along the trajectory. Typically, these 
are calculated with quantum chemical methods, such as 
time-dependent density functional theory (TDDFT) in 
Ref. [4] or semi-empirical methods in Refs. [BEIEIEIIITI. 
Mostly, these methods account for the influence of the 
protein on the pigments (here the BChls), by treating the 
protein environment as a collection of partial charges. 

An alternative to quantum chemical calculations are 
classical electrostatic calculations. The latter gave 
promising results for site energies, calculated using the 
crystal structure, without including dynamics m- It 
has been suggested in Ref. m to use this approach 
combined with MD simulations. Recently, this proposal 
has been implemented [3], calculating transition ener¬ 
gies using the so-called charge density coupling (CDC) 
method [18]. There, one calculates the transition energy 
from the interaction of the partial charges of the protein 
with ’transition charges’ of the BChls. 


As indicated above, the transition energy along a tra¬ 
jectory (i.e. the energy gap function) plays a fundamental 
role for both, the determination of the site energies and 
spectral densities. We will show in this article that the 
details of the molecular dynamics simulations (for ex¬ 
ample the force fields, or the initial condition) and the 
methods used to calculate the transition energies (clas¬ 
sical versus quantum chemical, and choice of quantum 
chemical method) strongly influence the obtained energy 
gap functions and the derived quantities. To make com¬ 
prehensive checks on convergence and initial state prepa¬ 
ration, we use computationally fast methods to calculate 
the electronic energies. Our results thus offer guidance 
in setting up the MD and energy calculation architecture 
for future applications of more advanced and computa¬ 
tionally intensive methods. 

The present work does not attempt a fully quantita¬ 
tive description of the FMO complex. Rather we wish to 
present an overview of the implications and limitations of 
the methods commonly employed for such a description. 
To facilitate this, we focus on the monomeric structure of 
the FMO complex. We include the 8’th BChl molecule 
despite it dissolving in MD simulations [T] , since the un¬ 
derlying structure determination [19] (PDB code 3BSD) 
is more recent than for the corresponding structures with 
7 BChls. 


The article is organized as follows: In section [TT] we 
will define the Hamiltonian underlying the open system 
framework for which we seek to extract parameters from 
the MD simulations. Then, in section HI details are given 
about the MD simulations, the calculations of the ’gap 
energies’ and determination of open system model pa¬ 
rameters from these. Section |IV| contains our results, 
grouped in basic properties of energy-gap trajectories 
dlv^ , e nergy correlations ( jlV C[ ) and pigment spectral 
densities (IVD). We summarize and interpret our results 
in the conclusion, section |V| A comprehensive collection 
of our data is provided in the supporting information [20]. 


II. OPEN SYSTEM APPROACH, SPECTRAL 
DENSITIES 


One often uses open system techniques to calculate 
optical- and transport properties of pigment protein com¬ 
plexes, as discussed in the introduction. To set the 
stage for these, consider a Hamiltonian divided into 
H = Hex + Hyii, + Hex-vih- Here Hex is the system 
Hamiltonian containing the relevant electronic states for 
the optical and transport properties. For linear spec¬ 
tra and excitation transfer it is sufficient to consider a 
restricted basis of states | n) = | 5^) • • • | e )n • • • | ^ ), in 
which BChl n is electronically excited to | e) while the 
other BChls are in their ground state \ g). Then we have 

Hex ~ n\n){n\+'^Vnm\n){m\ ( 1 ) 

n 


nm 
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where En is the transition energy, i.e. the energy dif¬ 
ference between ground and excited state of BChl n, 
and Vnm is the interaction associated with a transfer 
of excitation from BChl m to n. Throughout, indices 
n^m = 1... 8 label the BChl molecules as shown in 
Fig. [3b). 

The interaction has a strong dependence on the dis¬ 
tance and the orientation of the BChl molecules. In its 
simplest and most common form Vnm is written in the 
point dipole-dipole interaction m as 




^nm 


[fJjn ' l-^m ( 2 ) 


to determine thermal classical (electronic ground state) 
motion of the complex using molecular dynamics simula¬ 
tions described in section ITlIAl From the molecular coor¬ 
dinates thus obtained, we use semi-classical or quantum- 
chemical methods to obtain time-dependent energy dif¬ 
ferences between ground- and excited states of the BChls 
along the MD trajectory, see section [TlIB| From these, in 
turn, we extract the parameters of the model presented 
in section nn 


A. Molecular dynamics 


where fin is the transition dipole vector of BChl n and 
Rnm is the distance vector between the Mg atoms of 
BChls n and m. The factor / takes into account the di¬ 
electric environment, see iia[22] and references therein. 
We show Eq. for completeness, but do not consider 
interactions in this article. 

The environment is modelled as an (infinite) bath 
of harmonic oscillators HyUj = = 

E„ Ei hujnjaEan 

tion is taken as H. 


and the system-environment interac- 

x:-vih = Sn=l \^n){'^n \ (<^lj + 

Unj). Note that here we have assumed that all BChls 
couple to independent baths. This is probably a good 
approximation for internal vibrations of the BChls, while 
nuclear dynamics of the protein could simultaneously af¬ 
fect multiple BChls. In section [TV C| we will discuss this 
point further. 

The relevant quantity to describe the effect of the bath 
on one BChl molecule is the so-called spectral density. 


^ ^ l^njl ^(^ (3) 
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which describes how strong the electronic transition of 
BChl n couples to its oscillator j with frequency ujnj- 
Typically, these spectral densities are considered to be 
continuous functions of frequency. 

The SDs enter open quantum system approaches via 
the so-called bath-correlation functions (see e.g. Ref. O 

ESI) 


1 /> oo / \ 

Cn{t) = — duj Jn{uj) coth ( —- ) cosiujt) — i sm(ujt) 

TTJo L \‘^kBT J 

( 4 ) 

Here n refers again to the respective BChl and T is the 
temperature. 


All our calculations have been performed for one 
monomer of the FMO complex of C. tepidum for which 
we have taken the X-ray structure 3BSD from the Protein 
Data Bank as starting point. This structure is embed¬ 
ded in a water box with equal amounts of sodium and 
chloride ions added to obtain a 0.1 M/L salt concentra¬ 
tion. After energy minimization and equilibration, we 
have performed MD simulations of the ground state dy¬ 
namics using the MD simulation package NAMD m- 
We performed simulations using either the AMBER99 
[25] or the CHARMM27 [26l [27] force fields. Eor AM¬ 
BER the force field of the BChls was adopted from [28] 
and was also used in Refs. mn]. Eor CHARMM, we used 
the force field underlying Ref. [T] and has also been used 
recently in Ref. m- 

Eor each force field, we simulated molecular dynamics 
at different temperatures and with various initial condi¬ 
tions, described in detail below. Eor each initial condi¬ 
tion, we performed production runs of at least 300 ps 
duration, where the coordinates were saved every 5 fs. 
We describe further technical details of our MD simula¬ 
tions in m- 


1. Initial conditions for MD 

Here we briefly describe how we have generated differ¬ 
ent initial conditions (conformations) for the production 
run, in order to assess whether the usual assumption of 
ergodicity is justified in our context. To this end we used 
three different equilibration sequences (see [20] for more 
technical details on the equilibration step). 

A) We equilibrated for 10 ns at temperature T. 


III. METHODS 

Since the system part and the environmental part of 
the model introduced in the previous section both arise 
from molecular quantum dynamics of the EMO protein, 
they could in principle be extracted from ab-initio mod¬ 
elling of this structure. However, as a full quantum de¬ 
scription of the complex is out of the question, we resort 


B) We equilibrated for 10 ns at 300 K, then set the 
temperature to T and equilibrate there again for 
10 ns. 

C) We equilibrated for 10 ns at 300 K, heat up to 
310 K, then set the temperature to T and equi¬ 
librate there again for 10 ns. 

Eor T =300 K only sequence A was used. 
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B. Extraction of the open system parameters 

The output of each MD simulation described in the 
previous section is a trajectory containing the positions 
of all atoms in the FMO complex. 


1. Extraction of energy shifts 


In a first step, we calculate so-called gap-energies of the 
BChls, the energy difference between ground and excited 
state along this trajectory. They are written as 


En{t) — € + Sn{t), (5) 


where e is an arbitrary energy offset common to all BChls 
and Sn{t) are the deviations from this energy. From the 
En{t) we then calculate mean site energies Sn and spec¬ 
tral densities as detailed below. We will focus 

here entirely on differences between the gap-energies, for 
which calculation methods give reasonable results while 
absolute energies are problematic. Thus the offset e will 
drop out [36] . 

We compare two methods to calculate the shifts Sn{t): 

a) The first approach is based on classical electrostat¬ 
ics. It has been introduced in Ref. [18] for a station¬ 
ary situation and is denoted ’charge density cou¬ 
pling’ (CDC). In this approach, the energy shifts 
Sn{t) of the transition energy caused by the field 
of L point charges Qj in the surrounding dielectric 
medium are calculated as [18] 


(5 


n 


J_ ^ ^ Qj 


(6) 


Here Rj is the position of the jth atom of the pro¬ 
tein background, with charge Qj, and is the 
position of the ith atom of BChl n with transi¬ 
tion charge The latter rep¬ 

resents the change of charge density between the 
electronic ground- and excited state of the BChl, 
constructed as in [29]. The charges Qj of the pro¬ 
tein background are contained as parameters in the 
description of the respective MD force field, and all 
positions Rj are extracted from the MD tra¬ 
jectory. 

The effective dielectric constant eeff is used to de¬ 
scribe the screening and local field effects by the di¬ 
electric environment. In accordance with Ref. m 
we choose eeff = 2.5. We note that Ceff can also be 
treated as a fit parameter, determined by compar¬ 
ison with experimental results. 


spectroscopic properties (ZINDO/S) combined 
with the Configuration Interaction Singles method 
(ZINDO/S-CIS). For the calculation we used 
ORCA 2.9 [30]. The BChl molecule is treated 
quantum mechanical and the protein background 
within a sphere of 20 A radius around each BChl 
is taken into account via its point charges from 
the force field. The influence of the radius of the 
sphere is discussed in some detail in Ref. [ 7 ]. To 
further speed up the calculations the phytyl tail 
of the BChls was replaced by a CH3 group. For 
the calculation of the gap-energy the 10 highest 
occupied and the 10 lowest unoccupied orbitals are 
considered. This procedure is similar to the one 
used in Ref. [1]. 


2. Extraction of the electronic Hamiltonian and static 
disorder: 


In a second step, we convert time-dependent energy 
shifts Sn{t) thus obtained to the parameters that enter 
the open system model defined in section [Hj 

The site energies of the BChls are simply taken as 
the mean values of the energy-gap trajectories, i.e. 

£n = En{t) = e + 5n{t), (7) 


where the overbar denotes a time average over the full 
length of the MD trajectory. In a similar manner the ex¬ 
citation transfer elements Wm c ould be calculated from 
the MD trajectories by Vnm = ynm{E). In the present 
article, we concentrate on the calculation of the site en¬ 
ergies and resulting spectral densities, deferring a dis¬ 
cussion of these excitation transfer elements to a future 
publication. 

Besides the mean values, the trajectories En{t) also 
provide a distribution for the probability of a certain 
energy. Such distributions can be interpreted as static 


disorder m and are presented in section jlV A| 


3. Extraction of the spectral density: 

To determine a spectral density from the temporal fluc¬ 
tuations of energy shifts, we use the same formalism as in 
Ref. [5]. For BChl number n a classical gap-correlation 
function 

C^^{t,T) = {An{t + to)An{to)) , (8) 

is calculated, where (• • • )to denotes an average over 
to [37]. We indicate explicitly that the classical corre¬ 
lation function <© is obtained from a MD simulation at 
a specific temperature T. Further 


b) In the second approach we use the semi-empirical 
Zerner Intermediate Neglect of Differential 
Overlap (ZINDO) method with parameters for 


An{t) = Sn{t) - Sn{t) (9) 

is the energy gap fluctuation around the mean. 
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From this real correlation function the quantum spec¬ 
tral density is calculated via 

POO 

Jn{uj) = 2 ■ f{LO,T) dtC^\t,T)cos{L0t). (10) 

Jo 

To avoid artefacts in the numerical Fourier transform 
caused by the finite length of our correlation function, 
we multiply the classical gap correlation function by a 
Gaussian window function. In our case this Gaussian was 
chosen such that it leads to a convolution of the actual 
SD with a Gaussian of standard deviation 20 cm“^. 

Since the spectral density entering our open system 
model is a temperature independent function, the pref¬ 
actor /(cj, T) should ideally compensate any temperature 
dependence in Q. 

As has been theoretically shown in Ref. [5] , the choice 
/(co’,T) = 2 k bT appropriate for our harmonic envi¬ 
ronment in the open system approach and seems to be 
a reasonable choice in the case of the FMO complex. 
More details on the complex question how to construct 
a quantum spectral density (or equivalently a quantum 
bath correlation function) are discussed in Ref. [5]. 


IV. RESULTS 


Table H lists the different combinations of methods to 
perform molecular dynamics simulations and to extract 
excitonic Hamiltonians from the resulting data that we 
have compared here. The description of all results in this 
section refers to the method labelling in the table. 

In section [TV A| we first discuss general features of the 
site-energy trajectories resulting from all these methods. 
We then analyse the resulting mean site energies in sec¬ 
tion |Iv^ correlations of energies in section |IV G| and 


spectral densities in section IV D 


For each method the left panel contains a 5 ps segment 
of a trajectory, including all sampled data points at 5 fs 
intervals. They illustrate the fluctuations of the gap- 
energy on these short time scales. When comparing the 
CDG method (left column) with the ZINDO calculations 
(right column) it becomes evident that the fluctuations 
in the GDG method are roughly an order of magnitude 
smaller than those obtained with ZINDO. We will come 
back to this point several times in this work. Note that 
the absolute energy in the GDG method can be shifted 
arbitrarily, since it only delivers gap-energy differences, 
see section IIIIBII To provide a feeling for the slow long¬ 
time dynamics, we show in the middle panel for each 
method a 300 ps trajectory that was smoothed over win¬ 
dows of 5 ps [32]. The magnitude of the fluctuations in 
these curves is much smaller than that in the left pan¬ 
els. However, for the ZINDO it is again larger than for 
GDG. Finally, we show the distribution to have a certain 
gap-energy in the trajectory in the right panels, with a 
Gaussian fit as a guide to the eye. As expected from the 
left and middle panel, the width of the distributions for 
the ZINDO calculations is roughly an order of magnitude 
larger than that of the GDG calculations. 

We find that differences between these general features 
of the trajectories for different force fields (upper and 
lower row in Fig. are less pronounced than for differ¬ 
ent methods to calculate the gap-energies (left and right 
column). Note that MD simulations with different force 
fields employ different initial states and different trajec¬ 
tories. This apparently causes the ordering of the peaks 
in site-energy distributions of the GDG method (right 
panels, left column) to differ between the force fields. 
However, although we have used the same MD trajectory 
within each forcefield (i.e. within each row), the ordering 
of the peaks also differs between CDG and ZINDO. This 
comparison is made difficult by the large widths found 
with ZINDO, but is still possible e.g. for site 2. We will 
consider this point in more detail when discussing the 
site energies in section |rVB[ 


A. General features of the trajectories 

Before we discuss the parameters for the open system 
model it is worthwile to take a look at the trajectories ob¬ 
tained with the different approaches. They are shown in 
FIG.j^for all four methods at temperature 300K. Trajec¬ 
tories for the other BChls, at further temperatures and 
using different initial conditions can be found in the sup¬ 
porting information [20]. The trajectories presented in 
Fig. [^already allow us to discuss some general features. 


method 

la lb 2a 2b 

Forcefield 

Energy gap calc. 

CHARMM CHARMM AMBER AMBER 
CDC ZINDO/S CDC ZINDO/S 


TABLE I: Overview of calculation methods for molecular dy¬ 
namics simulations and BChl energy gap calculations used in 
the present work. 


1. Interpretation 

One can understand the difference in magnitude be¬ 
tween the GDC and the ZINDO calculations, when con¬ 
sidering the particular molecular properties that each 
method takes into account. In the ZINDO method, 
one calculates the energy difference between two Born- 
Oppenheimer surfaces. Here internal nuclear coordinates 
of the BGhls play an important role. Small changes away 
from the equilibrium position can lead to large changes 
in the gap-energy. This is sketched in FIG. [^, for the 
case of only one nuclear coordinate and shifted harmonic 
potentials. On the other hand nuclear coordinates ex¬ 
ternal to the BGhl affect the CDG only in very weak 
fashion. They only slighty modify the distances between 
the partial charges of the BChls and the protein. This is 
sketched in FIG. H). 
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BChl nr: — 1 — 2 — 3 


4 — 5 — 6 — 7 — 8 


13400 

13200 

13000 

12800 

12600 

12400 

12200 






[CHARMM/CDC| 




0 1 2 3 4 5 



3 100 

200 

[amber/cdc] 







12 3 4 

time in ps 


100 200 
time in ps 


100 200 
time in ps 


0 0.005 

probability 


FIG. 2: Example of the time-dependent site energies at temperature 300K. In each quadrant: (left panel) a zoom onto a 5ps 
segment, (middle panel) Full 300 ps trajectory smoothed over windows of 5 ps [32], which is the time interval shown in the left 
panel, (right panel) Probability to hnd a given gap energy within an interval of width 10 cm“^. The assignment of colors to 
pigment sites is shown in the legend and followed throughout this article. 


2. Lower temperatures 


It is instructive to repeat the calculations above for 
lower temperatures. We present results for 77K in Fig.|^ 
and results for 200K in the supplemental information 
[20] . The most striking new feature, is that the mag¬ 
nitude of gap-energy fluctuations, and thus the width of 
their distributions, is narrower than for the 300K trajec¬ 
tories. We find a roughly linear relationship of this width 
with the temperature T, for our three values of T. 

Similar to the 300K case, the energetic ordering of the 
different sites still varies between all methods for T = 11 
K, and is further different from the ordering at 300K. 
In contrast to room-temperature, we can clearly assign a 
relative order to the peaks of the energy-gap distributions 
also for the ZINDO calculations. Even when just varying 
the initial conditions as described in section |mA ^ we 
find variations in the peak ordering. 

While the width of the gap-energy distributions for 
all monomers are roughly of the same size at T = 300 
K, at T = 77 K some widths are nearly a factor two 
broader than others for ZINDO calculations. This holds 
in particular for site 6. We believe this is linked to the 
energy drift visible for site 6 in the middle panel of the 
upper right quadrant in Fig. at the beginning of the 
trajectory, and attribute it to conformation changes of 
the BChl, which cause a large change of the internal BChl 
energy-gap. Note that the CDC calculations in Fig. 
do not show the drift, although they are based on the 


same trajectory, presumably since CDC is less sensitive 
to internal BChl coordinates as argued in the previous 
section. 

Let us remark that Ref. [33], where a ZINDO approach 
(with CHARMM) has been used, also finds the site en¬ 
ergy width of some BChl to be much broader than of oth¬ 
ers. A trimeric complex has been studied in that work, 
which further finds that the energy width (and mean) of a 
certain BChl (e.g. no. 3) also changes drastically between 
the three monomers. 


B. Site energies 

From the trajectories presented in FIG. [2]|3] above, we 
now construct the mean transition energies of the BChls 
according to Eq. 0. The results are shown in Fig. § 
for CDC trajectories at 77 K, 200 K and 300 K, includ¬ 
ing variations for other initial conditions. We show the 
corresponding ZINDO results only in the supporting in¬ 
formation m, since the standard deviation of the mean 
transition energies there exceeds average energy differ¬ 
ences between the sites, complicating an interpretation. 

We can see in FIG. [^that all methods roughly agree 
on the overall site energy trend. The most striking dif¬ 
ference between CHARMM (method la, upper panel) 
and AMBER (method 2a, lower panel), is that sites 2 
and 5 have significantly lower energies using the former 
forcefield. Another important observation is that also 
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BChl nr: — 1 — 2 — 3 


4 — 5 — 6 — 7 — 8 




13400 - 
13200 - 

I 

13000 - 
12800 
12600 - 
12400 [ 
12200 - 
12000 ; 



13200 

13000 

12800 

12600 

12400 

12200 

12000 





2 3 

time in ps 


100 200 
time in ps 


0 0.05 

probability 



[CHARMM/ZINDO] 


100 200 




0 0.01 



12 3 4 

time in ps 


100 200 
time in ps 


0 0.01 
probability 


FIG. 3: Same as Fig. but for 77K. 



FIG. 4: We sketch how nuclear dynamics internal- or ex¬ 
ternal to the BGhl differently affect its gap-energy. (A) In¬ 
ternal vibrations. If the molecule is deformed the energy gap 
changes according to the difference of the ground and excited 
state Born-Oppenheimer surfaces. (B) External vibrations: 
The motion of the protein leads to changes in the electro¬ 
static interaction, which basically results in a energy shift of 
the difference between the two BO surfaces. In GDG only the 
latter effect is taken into account. Note, that protein dynam¬ 
ics can also induce dynamics within the BGhls. 


different initial conditions cause significant differences in 
mean energies. Since the underlying energy distributions 
are narrower than the differences between these means, 
as shown in FIG.|^ this indicates that the different equili¬ 
bration procedures that define our initial conditions (see 


section 


III A l| ) cause the system to end up in different 
We do not exclude that this ef- 


local energy minima, 
feet also occurs at 300K, but have not sampled different 
initial conditions there. 


Finally, we would like to point out that for some ini¬ 
tial conditions, neither site three nor four has the lowest 
gap-energy. These simulations then do not meet expecta¬ 
tions arising from the physiological function of the FMO, 
as transporting energy towards the reaction centre near 
sites three and four. However, most reported MD studies 
agree on site three having the lowest energy, as do most 
of our simulations. 


1. Comparison with previous results 


The literature values of the site energies of the vari¬ 
ous BChls show a large spread. Our results fall within 
this spread. Here we only discuss briefly the results of 
another CDC calcuation performed on the same species 
C.tepidum (but a different PDB structure) using the 
CHARMM forcefield [3]. Their values are quite similar 
to our AMBER/CDC results at 300K, the largest differ¬ 
ence being at site 6 where our value is roughly a factor 
1.8 larger. 

Our AMBER/CDC results also are close to the values 
of Ref. [15] with the notable difference that their BChl 1 
has a lower energy than BChl 2. 



































































FIG. 5: Mean site energies obtained from CDC simulations 
at 77K (triangles), 200K (squares) and 300K (circles) for the 
two different force fields. Upper panel: method la, lower 
panel: method 2a. The colors of symbols discriminate the 
different initial conditions (blue: A, red: B, orange: C). The 
color shading of the background identifies different sites as in 
FIG. El 


C. Correlations of site energy fluctuations 
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FIG. 6: Gorrelations Cnm according to Eq. (11) between the 


energy gap fluctuations at different sites n, m, using T = 
300iF. The diagonals fulfill Cnn = 1 and are set to zero. 
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In the definition of our open quantum system model in 
section [n] we have assumed independent harmonic baths 
for the different BChls. Although not necessary, this 
assumption makes calculations easier and is often used. 
Previous investigations indicate that indeed correlations 
between the gap-energy fluctuations at different BChls 
are very small [8]. 

In the following we investigate how strong the gap- 
energy fluctuations at different sites are correlated using 
all methods in table [H To this end we evaluate the cor¬ 
relation function between sites n and m 

^ _ ^n{t) ' ^m{t) /-j-jN 

^nm — : 

_ 1/9 

where C7„ = ((A„(i))2) is the standard deviation of 
the gap-energy distribution of BChl n, associated with 
the trajectory. Figs. and show the resulting correla¬ 
tions for 300K and 77K, respectively. 

Clearly, the correlations in the CDC method (left col- 


FIG. 7: Same as Fig. but for 77K. 

umn) are much larger than those obtained with the 
ZINDO calculations. Decreasing the temperature also 
decreases the correlations, an effect that is more pro¬ 
nounced for the CDC method. To quantify these obser¬ 
vations, we provide the maximal correlation, the sum of 
the correlations for all methods, temperatures and initial 
conditions in the supporting information [20] . 

The fact that the correlations for CDC are much 
stronger than for ZINDO is consistent with our obser¬ 
vations in section [IV A 1[ that ZINDO fully samples in¬ 
ternal potential energy surfaces, while they only have a 
weakened effect in CDC. Hence the main fluctuations in 
CDC are due to protein motion. Since it is fair to assume 
that protein motion can affect several BChl in a corre¬ 
lated fashion, while internal vibrations of separate BChls 
are uncorrelated, fluctuations that are more strongly af¬ 
fected by internal vibrations (ZINDO) are less correlated 
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between sites [38] . 

1. Comparision with previous results 

Our results here are consistent with a previous 
study on the FMO Chlorobaculuui tepiduui (PDB code 
3ENI) [8], which also essentially indicate the absence of 
correlations using ZINDO, with Cnm ^ 0.05 throughout. 
However, here we also see that correlations exist, if one 
is only interested in the external influence, that can be 
separated using CDC. 


SDs in CDC only indirectly via interactions with the pro¬ 
tein, in contrast to ZINDO. Since in both cases motion 
is determined by the same (ground state) energy surface, 
vibrational frequeneies of internal modes do however still 
affect the CDC results. Therefore we expect to see peaks 
in CDC spectra at similar positions as for ZINDO, but 
with much smaller magnitude. A clearly visible example 
is the peak near 1500 cm“^, present in both CDC and 
ZINDO. 


1. Speetral densities from CDC 


D. Spectral densities 


After our detailed discussion of features in gap-energy 
trajectories, such as the width of fluctuation distributions 
and the absence of inter-site correlations, we now proceed 
to the extraction of spectral densities 

We determine spectral densities from trajectories 
of energy-gap fluctuations according to Eq. and 
Eq. dTol) [39]. 


The spectral densities of all 8 Bchl molecules calculated 
with CDC (methods la and 2a) at 300K are shown in 
Fig. i Those calculated with ZINDO (methods lb and 
2b) at 300K are shown in Eig. In each panel we also 
provide integral measures that are often used to assess 
the ’total strength’ for the interaction of the respective 
BChl with the bath. The first measure is the so-called 
reorganization energy, defined as 


In general the SDs that we have obtained with CDC 
only show minor variations between different tempera¬ 
tures and different initial conditions. We emphazise this 
point, since for the ZINDO calculations we will see large 
variations in the magnitude. 

The results of EIG. can be interpreted as captur¬ 
ing the coupling to a bath that is external to the BChl, 
based on our discussion above. This has been more 
rigorously shown by Jing et al. [12], who applied CDC 
and ZINDO/S to MD simulations of the reaction cen¬ 
ter of purple bacteria. They compared ZINDO calcula¬ 
tions with and without partial charges on the surround¬ 
ing protein, and found that spectral densities based just 
on the differenee between these two fluctuations agree 
with those obtained from CDC. 


2. Speetral densities from ZINDO 


1 

Er = — duj J{(jj)/uj^ (12) 

^ J c 

the second one is the total Huang-Rhys factor 

x=- dL0J{w)/w'^. (13) 

^ J C 


We choose a small finite lower bound of c = 20cm ^ for 
the frequency integral, to avoid the divergence at cj = 0 
m- Here we exemplarily discuss the SDs obtained for 
300K, the SDs for other temperatures and varied initial 
conditions are shown in the supporting information [20] . 

Before we separately discuss detailed features of SDs 
obtained from CDC and ZINDO, we highlight the pro¬ 
nounced differences between the CDC and the ZINDO 
calculations, i.e. between Eig.[8]and Fig-H Two facts im¬ 
mediately catch the eye: Eirstly, the height of the peaks 
grows with increasing energy in the ZINDO calculation, 
while it roughly remains the same or even decreases in 
CDC. Secondly, the overall magnitude of SDs obtained in 
CDC calculations is much smaller than of those obtained 
using ZINDO. The latter, is also reflected in reorganisa¬ 
tion energies and total Huang-Rhys factors. This is again 
due to a distinction of internal vibrational modes of the 
BChls and external motion of the protein. As discussed 
in section IIV A 11 internal vibrations contribute to the 


The overall shape for SDs of all BChls is very simi¬ 
lar, independent of temperature and initial conditions. 
However, we notice two variations: 1) The highest en¬ 
ergy peak is sometimes strongly suppressed, for example 
in method 2b at 300K for site site three and four. 2) The 
magnitude of all SDs varies considerably from site to site 
as can be seen in Eig. [^ for 300K, and in the supple¬ 
mental information [20] for 77K and 200K. The smallest 
and largest SD differ by more than an order of magni¬ 
tude (as is also evident in the reorganization energy and 
the total Huang-Rhys factor). Eor a more quantitative 
picture, Eig. [To] shows reorganization energies for all tra¬ 
jectories analysed [41]. One sees large variations from 
site to site, within each temperature and for the different 
initial conditions. Nonetheless, as we discuss in the next 
section, the average of these results (per site) are roughly 
consistent with other theoretical studies and available ex¬ 
periments. The variations thus highlight that great care 
has to be taken in the detailed implementation of the 
schemes employed here. Our interpretation is as follows: 
The structure relaxation of the MD simulations yields a 
variety of BChl conformations with significant probabil¬ 
ity. The shift between ground- and excited energy sur¬ 
faces may depend on the conformation, yielding different 
amplitudes in spectral densities. Possibly more impor¬ 
tant, BChl conformations far from the equilibrium con- 
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FIG. 8: Spectral densities at 300K for method la (left) and method 2a (right), separately for each site. We also provide the 
values of the reorganization energy Er^ Eq. (12), and the integrated ’Huang-Rhys-factor’ X, Eq. (13). 


figuration of an isolated BChl might occur, for which the 
employed ZINDO method could give erroneous results. 
This may explain outliers in FIG.[^ as also conjectured 
in m- More clarity could be gained by a joint analysis 


of reorganisation energies and BChl geometries, which is 
beyond the scope of the present work. Another possible 
cause of the spread observed in FIG. 10 can be water 


molecules that enter the protein and reach the vicinity 
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FIG. 9: Same as FIG.[^ but using method lb (left) and method 2b (right). 


of BChls. Since we did not save the positions of water 
molecules in our simulations, we defer an investigation of 
this possibility to future work. 

Finally, let us point out that there is a pronounced 
shift between the high energy peak for the different force 


fields. For the CHARMM/ZINDO results this peak is 
located at an energy approximateley 200 cm“^ higher 
than for the AMBER/ZINDO. This is consitent with the 
results found in the previous literature. 
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FIG. 10: Reorganisation energies, Eq. |T^ , for all spectral 
densities calculated using the ZINDO method. Top: method 
lb, bottom: method 2b, symbols as in FIG. 
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3. Comparison to other theoretical calculations 


Let us first compare our results with theoretical stud¬ 
ies in the literature, which employed only a single of the 
methods compared here. We begin with Valleau et al. [5] , 
who present spectral densities that agree quite well with 
experiment and are independent of temperature. For a 
comparison, one however has to keep in mind that Ref. [5] 
models the slightly different structure of Prosthecochloris 
aestuarii (PDB ID; 3EOJ) g], containing only 7 BChl. 
Our results of methods lb and 2b are shown in Fig. 
together with those of Ref. [5] for monomer 3. We have 
averaged all SDs obtained from different temperatures 
and initial conditions. The resulting SDs are about a 
factor of two larger than reported in [5] . We see that the 
results of method 2b are more similar to [5], as expected 
since method 2b and Ref. [5] both rely on the same force 
field (AMBER). Also the width of the energy-gap distri¬ 
butions (see IQ]) is comparable. 

Another MD study of Prosthecochloris aestuarii was 
reported in [Tj, using CHARMM/ZINDO (method lb) 
Although they used a slightly different structure (PDB 



uj in cm 
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FIG. 11: Spectral densities of site 3, comparison of our re¬ 
sults and theory from [5]. Left: low frequency part. Right: 
full spectral density. The red and blue curves are an aver¬ 
age of all SDs that we have calculated for method lb and 
2b, respectively. The gray, dashed curve is the result of [5] 
multiplied by a factor of two. 


ID: 3ENI), their results for the monomer indicate similar 
trends as we observed here for the CHARMM/ZINDO 
calculations. In particular spectral densities have also 
a large overall magnitude and the highest energy peaks 
are roughly at the same position where we find them for 
method lb. 


V. SUMMARY AND CONCLUSIONS 

We have modelled the thermal motion of the FMO 
complex for C. tepidum using two different MD force- 
fields. From the resulting atomic trajectories, we have 
then extracted BChl gap-energies using again two dif¬ 
ferent methods, classical CDC, Ref. m, and quantum- 
chemical ZINDO/S-CIS within the ORCA package [30]. 
We can use the direct comparison of the four different 
sets of result to elucidate the origin of the rather large 
spread of literature results HHsi regarding the mean val¬ 
ues of these BChl gap-energies, their relative ordering 
and the associated spectral densities. 

Regarding the distribution of energies, we find that 
CDC gives very narrow distributions (with width of the 
order of a few 10th of wavenumbers) compared to ZINDO 
which gives quite wide distributions (width of the order 
of several hundred wavenumbers). This can be under¬ 
stood, since a change in the bond-length within a BChl 
leads typically to much larger shifts of transition ener¬ 
gies than changes in the external protein positions. In 
the quantum chemical method ZINDO, the former effect 
is fully accounted for, but not in CDC. 

We confirm earlier results of Ref. [8] reporting very 
weak inter-site energy correlations when using ZINDO, 
but find that correlations are much larger when using 
CDC. This is again consistent with a picture that en¬ 
ergy fluctuations in CDC are dominantly determined by 
motion of the protein external to the BChls. These may 
influence several BChl due to the long range nature of in¬ 
teractions (and the neglect of screening). In the quantum 
chemical ZINDO calculations these correlations should 
still be present, but are there overwhelmed by the larger 
fluctuations due to seemingly uncorrelated internal BChl 
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vibrations. 

Also the spectral densities of CDC and ZINDO dif¬ 
fer substantially, those from CDC having a much smaller 
amplitude and missing high frequency peaks, consistent 
with our interpretation that these peaks come from inter¬ 
nal vibrations. Only the ZINDO spectral densities have 
amplitudes of the same order as obtained from experi¬ 
mental FLN results m- 

We also found pronounced differences between the two 
forcefields. Notably, we have obtained quite different site 
energies, and for the ZINDO calculations the main peak 
of the SD was at significantly different frequencies. The 
latter finding indicates that the forcefields that we have 
used for the BChls result in quite different high energy 
vibrations. 

When we averaged the SDs for a given force-field over 
the various temperatures and initial conditions of our 
ZINDO calculations, we found good agreement with ex¬ 
periments [14] and theory [5|, indicating rough consis¬ 
tency of the entire calculation scheme. However it is the 
large spread in basic outcomes of the calculation without 
this averaging, exemplified in FIG. 10, that we consider 
one central result of the present work. It indicates that 
great care has to be taken in all details of implement¬ 
ing the calculation scheme. Our interpretation of these 
variations, is that multiple local energy minima play a 


role when the BChl relax during MD structure equilibra¬ 
tion. In particular at low temperatures the system may 
stay for long times in a given minimum, but move be¬ 
tween different configurations on even longer time scales. 
Jumps between certain energy values have been observed 
in spectral diffusion experiments on single light harvest¬ 
ing complexes of various species [34l [35| . 

It seems that these different configurations have lit¬ 
tle influence on the SDs in the CDC approach but mas¬ 
sively affect those in ZINDO. An additional problem may 
be that some BChl configurations that are relevant for 
the MD are sufficiently far from the isolated equilibrium 
structure to cause problems for ZINDO. All these obser¬ 
vations necessitate a very careful MD equilibration and 
choice of quantum chemical method, and point to the 
need for further studies for a more complete picture of 
the interplay of these methods. 
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ternal vibrations of different BChls are uncorrelated and 
that they are also uncorr ela ted with Protein -dynamics. 
Then we find An(t)Am(t) ~ AB’^®^(t) AB’^°^(t) and finally 


Cnm ~ ^ jn \ This is small as one sees 

VABChl(^)2^^BChl(^)2 

using similar aguments as above. For CDC, assuming 
one has in 

the denominator of that expression. Thus Cnm will be 
larger as for the ZINDO case. 

[39] Convergence regarding the time-average in Eq. ^ was 
confirmed by comparing results obtained from the full 
300 ps trajectories with those from averaging over ran¬ 
domly chosen 5ps sub-segments, finding identical results. 

[40] The choice of c strongly affects X and Er. We choose 
it as the lowest value for which, increasing c does not 
significantly alter X. Note also that the correct descrip¬ 
tion of frequencies near zero would require much longer 
trajectories, such that ignoring these frequencies is more 
physical. 

[41] We obtain the same qualitative picture when plotting the 
maximal peak height or the total Huang-Rhys-factor. 

















